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ABSTRACT 



A computer model is constructed for the purpose of 
simulating interception of a near earth satellite by a 
grooind launched ballistic interceptor with emphasis on 
adaptability of the model to the investigation of interceptor 
performance in terms of interceptor and satellite character- 
istics. The model is keyed to instantaneous satellite 
positions, and provides on call information about minimum 
interceptor velocity requirements, parameters of possible 
interceptor trajectories and relative velocities at inter- 
ception for those trajectories. The model is constructed 
with sufficient generality and flexibility so that it may be 
readily adapted to simulation or computation of other 
intercept problems. 
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CHAPTER I 



INTRODUCTION 

The fundamental equations of motion for a body 
acting under the influence of a central force have been 
•well known and exhaustively treated since the times of 
Newton and Kepler. These equations, which will not be 
redferived here, form the basis of the model to be devel- 
oped. It should be noted, however, that this simple 
formulation involves some assumptions that do not com- 
pletely correspond to the physical environment of a near 
earth satellite. These assumptions and their effects on 
the model must be considered. 

BASIC ASSUMPTIONS 

1. The earth and the satellite are considered as point 
masses with the mass of the satellite negligible compared 
to the mass of the earth. 

2. The earth's gravitation is the only force acting on 
the satellite. 

3 . The coordinate system is fixed in space. 

DISCUSSION OP ASSUMPTIONS 

If the earth were perfectly spherical with its mass 
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distributed symmetrically^ about its center 5 all 'its mass 
could be considered to be concentrated at one point. 
However the earth is slightly oblate to a first approx- 
imation and there are slight asymmetries in the distri- 
bution of its mass, particularly in the surface distribution 
of continental and ocean masses. This oblateness is 
responsible for the major anomalies in its gravitational 
field, which affect the behavior of near earth satellites 
in two important ways [3] . First, it causes the orbital 
plane to rotate about the earth's axis in the direction 
opposite to the satellite's motion. This perturbation is 
known as regression of the nodes and its rate is given by 



- 9.97 (— 

Qu 



3.5 



2 —2 

( l-e ) cos i degrees/day (1.1) 



where R is the radius of the earth, a is the semi-major 
diameter of the orbit, e is the eccentricity of the orbit 
and i is the inclination of the orbital plane to the plane 
of the earth's equator. These terms and symbols are 
discussed in detail in Chapter II. The notation through- 
out this paper will be in accordance with the table of 
symbols vinless otherwise specified. Second, the earth's 
oblateness causes the perigee position of the orbit to 
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rotate about the principal focus of the orbit with the 
rate being given by 



d6? ~ If. 98 



3 3 3 3 3 

(l-'O ) (5 cos i-1) degrees/day. (1.2) 



dt la 

The atmosphere offers increasing resistance to the 
motion of a satellite as its altitude decreases. The 
primary effect of atmospheric drag is to reduce the major 
diameter of the orbit, mostly at the expense of apogee 
distance. This makes the orbit more nearly circular with 
the passage of time while reducing the satellite's period 
of revolution and increasing its velocity in accordance with 
Kepler's third law. 

In addition to these major effects, satellite orbits 
are perturbed to a much smaller degree by solar wind, 
radiation pressiire and the gravitational attraction of 
other celestial bodies, especially the moon. 

For most applications of the model to be developed, 
the errors mentioned above may be disregarded. The 
model is geared to inspection of instantaneous satellite 
position and velocity vectors, and orbital parameters may 
be considered to be no older than one period of revolution. 
However, the model should not be applied to orbits having 
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altitudes of less than 100 miles because of the increased 



effects of atmospheric drag, nor to orbits with apogee 
distances of several earth radii because of perturbations 
due to the moon's gravitational influence. Neither should 
the model be used to predict the position of the satellite 
for more than one or two orbital periods because of the 
cumulative errors indicated by equations (1.1) and (1.2). 
These restrictions are not severe, for in fact most 
objects of interest in the interception problem, such as 
reconnaissance and weapon carrying satellites will orbit at 
altitudes between 100 and 1000 miles due to the require- 
ments of their missions. Substitution of these values 
into equations (1.1) and (1*2) reveal that even in the 
worst cases, the perturbations amount to only a fraction 
of a degree per revoluj:ion. 
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CHAPTER II 



ORBITAL MECHANICS 

Kepler's brilliant formulation of the laws of planetary- 
motion early in the seventeenth century were provided 
with a sound theoretical basis by Newton's laws of motion 
and were shown to apply generally to bodies acting under 
the influence of a central force. These laws can be 
stated for an earth satellite as follows: 

1. The orbit of a satellite is an ellipse with the 
center of the earth at one of its foci, 

2. The radius vector of a satellite describes equal 
areas in equal times. 

3 . The square of the period of a satellite is 
proportional to the cube of its mean distance from the 
center of the earth, 

THE ELLIPTICAL ORBIT 

The size and shape of the elliptical orbit is deter- 
mined completely by the parameters a and e. (See Pig. 
2-1). These quantities can be expressed in terms of 
other orbital dimensions when a and e are not given 
directly. For example, the ellipse is determined if only 
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ap ogee 




perig ee 



a semi-major diameter 

b semi- minor diameter 

C center of the earth 

e eccentricity 

r satellite position vector 
p perigee distance 

argument of the perigee 
Q true anomaly 

X direction of ascending node and x-axis 

L semi-latus rectum 



Figure 2-1, The Elliptical Orbit 
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the minimum and maximum altitudes of the satellite are 



known. The perigee distance p is given by the minimum 

altitude h . plus the radius of the earth R, Similarly, 
min ' 

the apogee distance q is given by h + R. The semi- 
max 

major diameter a is then (p+q)/2 and the eccentricity e, 
the semi-minor diameter b and semi-latus rectum L are 



given respectively by 



■\Ti 2 

= a V l_e 






( 2 - 1 ) 

( 2 - 2 ) 

(2-3) 



a 

from the geometric properties of the ellipse. The angular 
displacement S of the satellite from the perigee is 
measxired in the direction of satellite motion and is 
commonly called the true anomaly. 

It remains to determine the position of the satellite 
in its orbit as a fiinction of time. In polar coordinates, 
this position is given by 

r a(l-e^) (2-4) 

1+e cos*0 

Unforttmately 0 is not a simple function of time, 
and the satellite position is more easily expressed in 
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terms of the eccentric anomaly E, whose geometric 
relationship to the satellite position is illustrated in Fig. 

2 - 2 . 



From Newton's laws, the radial acceleration 
2 

(9»80665 m/sec ) of a body in circiilar orbit of radius R 
is equal to /R. Solving for v^ yields a velocity of 
7,904*455 meters/second. 

Then the period in seconds is given by 

P = 27/^ = 5104.6536 seconds. (2-5) 

O V 

''o 

By Kepler's third law, the period of any other near earth 
satellite is then 



a\ 3/2 

P = P '—I seconds, 
and the mean angular motion p. is given by 
p. = radians/second. 

Now define the mean anomaly as follows. 



(2-6) 



(2-7) 



M = ^ (t-T) 



( 2 - 8 ) 



where t is time in seconds and T is the time the satellite 
passes throTigh perigee. From this, Kepler's eqiiation 

M = E -e sin E (2-9) 

can be solved for E, (See Chapter III). From Fig. 2-2, 
it is apparent that 

r == a (1 -e cos E) (2-10) 
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tulo* folo* 



y' 




= (1-e cos E)^ 



Figure 2-2. Position vector in terms of E 
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from which, together with equation {2-4) 



cos 6 



cos E -e 
i -e cos E 



sin 9 = -Vl -e^ sin E 






( 2 - 11 ) 



Equations (2-11) make it possible to write the cartesian 
coordinates of S in Fig. 2-2 as 



x' = r cos 9 = a (Cos E -e) 
y’ = r sin 6 = a Vl -e^ sin E, 

j 

GEOCENTRIC EQUATORIAL COORDINATES 



( 2 - 12 ) 



The reference system to which the orbital parameters 
are referred is known as the geocentric equatorial 
coordinate system. It is a conventional right-handed 
cartesian system with the origin at the center of the 
earth. (See Fig. 2-3). The x-axis is directed toward 
a fixed point in space where the sun's apparent path 
through the heavens crosses the earth’s equatorial plane 
in a northerly direction. This point is variously known 
as the first point of Aries or the position of the vernal 
eqmnox. It is not truly fixed in space as stated in the 
basic assumptions, but its movement due to nutation of 
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z 




nr direction of vernal equinox 

jTL longitude of ascending node 

CJ argument of perigee 

i inclination of orbit 

0 true anomaly 



Figure 2-3. Projection of Orbit on Unit 
Geocentric Sphere 
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the earth's spin axis is only a small fraction of a degree 
per year and can be ignored in this problem. The y-axis 
is 90° from the x-axis, measured eastward along the 
equator, and the s-axis is in the direction of the north 
pole. Corresponding unit coordinate vectors are designated 

T, J and k. 

The orientation of the orbit plane in the reference 
system is given by the parameters jO , cJ , and i. The 
longitude of the ascending node fl is the angle measxired 
eastward along the equator from the x-axis to the point 
where the satellite crosses the equator plane in a northerly 
direction. The argument of the perigee O is the angle 
measvired in the orbit plane in the direction of motion 
from the ascending node to the perigee position. The 
inclination i is the angular separation of the equator and 
orbit planes. This angle is less than 90® if the satellite 
motion is eastward (or direct). If the motion is west- 
ward (retrograde) the angle is the supplement of the 
angular separation of the planes. 

The change from x', y' coordinates in the orbital 
plane to x,y,z coordinates can be accomplished by the 
linear transformation 
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(2-13) 



X = 


^x 


X* 


+ 


Qx 


y’ 


y = 




X* 


+ 


^y 


y' 


z = 


^z 


x' 


+ 


Qz 


y’ 



\ 






/ 



where the P's and Q's are the direction cosines of the x' 
and y' axes relative to the x-axis [sj , These six 
quantities are found to be 



P 

Q 



X 

X 




P 

Q 



z 

z 



cos 63 cos Jl. - sin CJ 
-sin cO COS-0- — cos CJ 
cos cJ sinJl + sin 6J 
-sin cJ sinXL + cos UJ 



sin .a COS i 
sin D. COS i 
cos cos i 
cos JL cos i 



sin iJ sin i 
cos u) sin i. 



\ 



(2-14) 



/ 



The eccentric anomaly may be introduced directly 
into equations ( 2 - 13 ) • If we define 



Ax * 


B^ = aVi 




\ 






Ay. = aPy, 


By. = aV 1 -e^ 








(2-15) 


^z = 


Bg = aVl -e^ 


^Z 


> 







then 

X = (cos E -e) + sin E 

y = A ( cos E -e) + B sin E 
j y 

z = A_ (cos E -e) + B sin E 

z z 



( 2 - 16 ) 
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The first time derivatives of these position 
coordinates yield the components of the satellite velocity 
vector in terms of E. They are found to be 



dx _ 
dt 


r 


(-Ax 


sin 


E 


+ 


^x 


COS 


E) 1 


1 




^ = 
dt 


r 


(-Ay 


sin 


E 


+ 




cos 


E) 




(2-17) 


dz ^ 

■ar 


r 




sin 


E 


+ 


^z 


cos 


E) ^ 







CONSERVATION OF ENERGY. MOMENTUM AND 
AREAL VELOCITY 



Kepler's second law and the laws of conservation of 
energy and angiilar momentum reveal relationships which 
will be found useful. The famoiis Vis Viva integral, the 
energy equation for an elliptical orbit, states that 
/ \ 



v^ = GM 



r 



a 



P 2 

meters'^/ second 



(2-18) 



where G is a gravitational constant and M is the mass of 
the earth. In the MKS system of measurement, the 
product GM is equal to 3»98 x meters-^/seconds* and 

if nautical miles are used as units of distance, GM is 
eqioal to 6.26? x 10^ n.m.'^/seconds'^. It is noticed 
that for any orbit with a given major axis, v is a function 
only of r, or eqx 3 lvalantly,_ of altitjade. 
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Angular momentum is conserved in the two body 



problem, and its vector representation is given by the 
constant vector 



H = r X 



dr 

dt 



(2-19) 



Multiplication by r. yields the equation of the orbital plane. 



r • H = O. 



( 2 - 20 ) 



From Kepler's second law, the areal velocity of a body in 
an elliptical orbit is constant [ 2 ] . It can be shown that 

= -L H. (2-21) 

dt 2 

Then by taking the time integral of both sides, the area 
swept out in time t is given by 

A= -i- Ht + constant. (2-22) 

2 

In particular, if the area swept out between two points r^ 
and r^ is known, than the time in the orbit between these 
two points is 

^21 u 21 - 
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CHAPTER III 



SATELLITE ORBIT COMPUTATION 

SELECTION OF ORBITAL ELEMENTS 

A satellite orbit may be determined from a wide 
variety of given parameters. In the present model, 
those initial conditions and parameters are set which the 
analyst may consider to be of particular importance. Of 
first interest are the maximum and minimum altitudes, 
from which the size and shape of the orbit are determined. 
Specifying the inclination immediately determines the 
northern and southern limits of latitude over which the 
satellite may pass and the direction of motion. Setting 
the argument of the perigee fixes the perigee point in 
the desired hemisphere, north or south. Finally, the 
time at perigee and longitude of the ascending node are 
chosen or computed as indicated below and the orbit is 
computed from the equations of Chapter I, 

POSITION VECTORS IN TERMS OF TERRESTIAL 
LATITUDE AND LONGITUDE 

The analyst may wish to know the components of a 
position vector corresponding to a position on earth 
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designated by latitude and longitude at a given time. For 
example, he may wish to specify an initial geographic 
position for the perigee of an orbit or determine the 
position vector of the interceptor launch site. Conversely, 
he may wish to know the geographic location of a satellite 
position when its position is given as a vector. 

The Greenwich Hour Angle of the first point of Aries 
(GHAT* ) is listed in the Nautical Almanac for the begin- 
ning of each day of the year [l] . GHA'IT is defined as 
the central angle measured westward along the equator 
from the Greenwich meridian to the vernal equinox. 

Designate the longitude of a position by \ with eastern 
longitudes being positive, and the latitude by S with 
northern latitudes being positive. Then the Sidereal Hour 
Angle (SHA) between the vernal eqviinox and the position's 
meridian is equal to GHaT +A. However SHA increases 
at the rate of the earth's rotation, which in the chosen 
coordinate system is equal to one revolution per sidereal 
day or .0000729 radians/second. If the value of GHaIQ 
is taken as that given for O January 1 of the year, and 
t is the number of seconds after the beginning of the year, 
SHA = GHA*Y^ + A + 0.000729t radians (3-D 
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and 



X = r cos SHA cos S 

y = r sin SHA cos S (3-2) 

2 = r sin S ^ 

If the position vector r(x,y, 2 ) is given the process 
may be reversed although care must be taken in assigning 
signs to the computed angles. 

(3-3) 
(3-4) 



,222 
'x +y +2 



r - a/ 

S = Sin 
SHA = sin"^ 



-Mf 



ir 



rcoSjj^l (3-5) 

where - ^ ^ ^ with the same sign as z and 

-rf ^ SHA ^ tt with its magnitude and sign determined 



by seeing which quadrant contains the x,y coordinates of 
r. \ may then be found directly from equation (3-1)* 

If, as suggested above, the analyst wishes to specify 
an initial perigee point ajid longitude of ascending node from 
geographical coordinates, the corresponding position vectors 
r^ and are determined from equations (3-1) and (3-2). 
XL is of course simply equal to SHA and 



CJ = cos 



-1 



- \ 

,-rt Y» J 



i 

T,. r , ; 






(3-6) 

with the same sign as 2 . It must be noted that these 
two points fix the orbital plane so that the inclination may 
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no longer be chosen arbitrarily. However, it is easily 



computed, for the vector 






(3-7) 



is inclined to the z-axis by the angle i if the satellite 

IT' 

motion is direct and — + i if the motion is retrograde. 



Then 



1 (or -^+ i) = cos |-~ 1* 



r; 



SOLUTION OF KEPLER'S EQUATION 



Given the orbital parameters and the satellite's 
position in the orbital plane, the computation of coordinates 
in the geocentric equatorial coordinate system is a straight- 
forward process. However, the position of the satellite 
in its orbital plane as a function of time is found from 
the solution of Kepler's equation which is transcendental 
in E. This equation is solved by an iterative method of 
successive approximations. 

Write equation (2-9) as 

E = M + e sin E^ (3-8) 

and setting E^ equal to M, solve for E. Next setting E^ 
eqtial to the computed value of E, repeat the computation. 
The procedure is repeated until the desired accuracy is 



19 



obtained. It can be shown that the iterated values of 



E converge rapidly to a single value. 
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CHAPTER IV 



INTERCEPTOR TRAJECTORY COMPUTATION 

The first problem to be considered may be stated as 
follows. If an interceptor laimch position and a predicted 
intercept point are given, how can one determine trajectory 
parameters that will yield a successfixl interception. In 
particular, out of all possible interceptor trajectories, 
which one requires the least energy to be supplied to the 
interceptor vehicle. 

After the interceptor has exhausted its propulsive 
energy, it is subject only to the same central force of 
gravitation as is the satellite. It will follow an -elliptical 
orbit with a major semi-diameter determined by its altitude 
and velocity in accordance with eqviation (2-18). Conversely, 
if the semi-diameter of the interceptor orbit is known, ^ 
the energy of the orbit is fixed. Thus the orbit with 
the minimum semi-diameter which passes through the launch 
site and the intercept position is the minimum energy 
interceptor trajectory. 

This reasoning assumes that the interceptor attains 
its maximum velocity immediately after its launch and 
thereafter follows a purely ballistic path with air resistance 
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neglected. Actually there is an initial phase of its flight 
during which it acquires the velocity and attitude 
necessary for the proper ballistic path to the intercept 
position. This phase depends upon the thrust character- 
istics of the propulsion system, the payload weight and 
other factors of interest to the designer and engineer. 
This paper considers only the ballistic reqviirements which 
the propiilsion system must meet, and the interceptor is 
taken to be injected into its orbit with the necessary 
velocity vector at the time of latinch. 

Consider the plane determined by the position vector 
of the launch site and the position vector of the intercept 
point as shown in figure 

Let C be the midpoint of line BI and let the 
distance CF= 1/2 GI . (Notation: BI represents the 
magnitude of the directed line segment BI.) Note that 
by construction the distance OG+GI+IF = OB+BF» Now 
define e^^ as the ratio Cl/ CF and Pj^ as the quantity 
CF(ej^^ -1). Also let bj^^ = CF(pj^). Construct the 
hyperbola HH' from the equation 

"•h =_Ph 

1 - ej^ cos 0 



22 



I 



» 






I intercept point 

R radius of the earth 

h satellite altitude 



Figure 4-1 • Locus of Secondary Foci of 
Interceptor Trajectories 
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with the origin at I. B is the focus of a mirror hyper- 
bola and from the geometric properties of the hyperbola, 
the distances from I and B to any point on the cxxrve 
( F' for example) differ by a constant, in this case eqiial 
to GIo Thus the distance OB+BF' is equal to the 
distance OI+IF'. But recall that the sum of the distances 
from any point on an ellipse to the two foci is a constant, 
equal to the major diameter 2a. Therefore any allipse 
with one focus at the center of the earth and the other 
lying on the hyperbola HH' will pass through the launch 
site and the point I . Clearly the ellipse with the smallest 
value of a is the one where F=F', thus determining the 
minimum energy orbit. 

The orbital elements of the minimum energy inter- 
ceptor trajectory may now be fo\md from the given 
position vectors OB and OI . BI=OI - OB. The distance 
GI is the altitude h of the satellite at intercept and is 



to OI — 


R. 




BF = i 


fBI+h_ \ gj 


(4-2) 


BI 


OF = OB + BF 


(4-3) 


a = OB + BF 


(4-4) 



24 



e 



OF 

2a 



( 4 - 5 ) 



The vector product OK = OB x OI is normal to the plane 
of the orbit, and its component will be positive if the 

motion of the interceptor is direct. Then 



1 = cos 



-1 



OK , k 
OK 



( 4 - 6 ) 



The vector product ON = k x OK points in the direction 
of the ascending node, so 



il. = cos 



-1 



ON . i 



\ 



ON 



( 4 - 7 ) 



Perigee distance is equal to a - OF, so the perigee position 
vector 



OP = (OF -a) OF 



( 4 - 8 ) 



and the final spatial element 



0 > 



-1 

cos 



ON . OP 

(on)(op) 



( 4 - 9 ) 



The period of the orbit and the mean motion p. of 
the interceptor are found from equations (2-6) and (2-7) » 
The sine of the true anomaly of the interceptor at the 
launch site is given by 
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sin 



e = 



OP X OB I 
OP) (OB) 



( 4 - 10 ) 



and from equation (2-11), sin E and E are determined. 
Equation (2-9) gives the value of M and finally 



t -T =~ (4-11) 

M 

from which the epoch of the orbit is determined. A 
similar computation will give the time from perigee to 
intercept from which the time of flight of the interceptor 
is known. 

The above derivation determines the orbital parameters 
for the minimum energy orbit but the same procedure may 
be followed whenever the location of the secondary focus 
is given. The orbital plane is determined by the location 
of the launch site and the intercept point, and any secondary 
focus position in that plane can be computed by equation 
(4-1). A transformation of coordinates will yield the 
position vector of F* from which the desired parameters 
may be computed by the methods above. 

The velocity vector can now be computed by equation 
(2-17) for any point in the orbit. Of couj?se, the points 
of greatest interest are the point of interception and the 
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0^ ■• villi: 









launch site. At the intercept point, the difference 
between the satellite and interceptor velocity vectors is 
the closing velocity. The velocity vector at the launch 
site is the sum of the velocities imparted by the propulsion 
system and the velocity vector of the launch site due to 
the earth’s rotation. At lower latitudes, this latter 
velocity approaches 1000 mph and an interceptor with a 
given propulsion system will have a substantially longer 
maximum range when launched to the east than when 
launched toward the west. The velocity vector of the 
launch site is seen from equation (3-2) to be 



dx 

dt 

dy 

dt 



-uR sin SHA cos S 
uR cos SHA cos ^ 



dz _ 
dt “ 



0 



( 4 - 12 ) 



This quantity must be subtracted from the orbital velocity 
vector at that point if the velocity to be supplied by the 
propulsion system alone is to be determined. 

Figure 4-2 shows the minimum energy orbit for the 
conditions displayed in Figure 4-1 together with another 
interception orbit having a secondary focus F' on the 
hyperbola HH'. 
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Figure J+-2, The Minimum Energy Trajectory 
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CHAPTER V 



COMPUTER PROGRAMMING CONSIDERATIONS 

The satellite interception problem, with its complex 
and constantly changing relationships of the many para- 
meters and variables involved lends itself very well to 
simulation by a modern high speed computer. Huge amounts 
of information can be computed, stored and recalled with 
great speed and accuracy. Such a computer model has 
been constructed for execution on a Control Data Corp- 
oration computer model IbOlf* (See Appendix A), 

Substantial economies can be effected by writing 
computer programs in machine language in terms of 
execution time and memory requirements. This comes at 
the expense of using a format which is not easily adapted 
for vise by other computers and is not readily modified by 
other programmers. An alternative, which is the approach 
used in this model, is to write the computer program in a 
commonly used compiler language, a variation of FORTRAN 
known as FORTRAN-63 • The model follows the 
mathematical approach outlined in this paper, and by 
noting the liberal number of explanatory comments in the 
program, an analyst should be able to adapt the program 
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for his own purposes on another computer employing a 
FORTRAN compiler langviage. The remainder of this 
chapter is devoted to inspection of computer and compiler 
langxiage considerations which affect the construction of 
the model. 

Inasmuch as the satellite orbit is fixed in the 
chosen coordinate system, much computational time can 
be saved by computing position and velocity vectors for 
one orbital period and storing the information for recall 
as needed. If sufficient storage space is available in the 
computer memory, all quantities of interest such as 
position vector components, velocity vector components, 
true and eqcentric anomalies with their trigonometric 
fxinctions, distance from the center of the earth, etc., 
may be stored for later use, and this is recommended 
where possible. However, the period of a circiilar orbit 
at 1000 miles altitude is nearly 7500 seconds, and storage 
of all quantities of interest in sufficiently small increments 
of time imposes a considerable requirement for memory 
space. 

In this particular model, only the position vector 
components can be stored for each second of time in 
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normal format -without exhausting the memory storage 
capacity. By going outside the FORTRAN -63 language 
to a compatible assembly program CODAP-1, it is 
possible to store the eccentric anomaly for each second 
by storing two quantities in a single storage address • 
The CODAP— 1 subroutine BYTE5 which accomplishes this 
is not included in the FORTRAN -63 listing, 

A small bank of temporary variable names (Cl through 
C9) in a COMMON statement saves some storage by using 
the same names wherever an intermediate, non-permanent 
quantity is computed. The same statement contains 
conversion qmntities that are needed in several sub- 
routines, Additional information may be stored on magnetic 
tape or equivalant devices but this is profitable only if the 
time required for access to the information is less than 
the time required to compute it. 

Liberal use is made of subrou-tines and functions in 
•the interest of flexibility, for the model is constructed 
so that more relationships than are explicitly specified may 
be considered with only minor modifications. For example, 
the model can easily consider a terminal guidance aspect 
of the problem or serve as the basis of a probabalistic 
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approach to interception. Communication between the 
subroutines is facilitated by the use of COMMON state- 
ments. 

Clarity of the program for other analysts has been 
enhanced by consistent use of sixffixes and prefixes in 
variable names and the use of mnemonic variable names 
throughout. 
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APPENDIX A 



fCrtran-63 computer. sijc^ulation program 

The program listed here has been extensively tested 
for a wide variety of orbital elements and the generated 
orbits show close agreement with the orbits of actual 
satellites. Storage is provided for satellite positions in one 
second increments for orbital altitudes to 1000 miles. 

Sample output statements have been included to suggest 
ways in which the user might call for computed values. 

The computational procedure follows the methods 
outlined in the thesis with extensive use of subroutines and 
functions for those calculations that are used repeatedly 
or in other sections of the program . This enhances the 
flexibility of the program for the analyst who may want 
the calculated quantities for supplementary considerations 
in the interception problem. 
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FUNCTION RHOF(X,Y,Z) 

FUNCTION VALUE IS THE MODULUS OF 
RHO = SCRTF(X*X + Y*Y + Z»Z) 
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412H (DEGREES) ,6X,gH( DEGREES) , 6X , 9H ( DEGR E ES ) ,7X,7H(KILES),8X,7h(M 
5 IL ES) , 19X , )6H{DAY/HR/MIN/SEC) ,2X,9H(M IN/SEC ),/ ) 

I = SATPER/60 $ J = XMODF( SATPER,60 ) 

PRINT 1080, SATNOD , S AT ARG , SA T I NC , SATHK X , SA THPN , S ATECC , I CAY , I HR, 
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COMMON /SATCCOR/ SAT XO, SATX (71*36) , SAT YO , SATY ( 7*4 36 ) , SAT20 , SATZ ( 74 36 ) 
COMMQN/CON/PI ,TV*OPI ,RAD,C1 ,C2,C3,C4,C5,C6,C7,C8,C9 
COMMON/ SA TP AR AM /S A TARGfSATNOC, SATING, SATHMN fSATHMX , SAT PER, SAT ECC, 
1 SATSMA, SATSMB ,AX, A Y , A 2 , B X , B Y , B Z 
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